############################################
# Author   : Aastha Rajan 
# Date     : 04/07/2020 
# Project  : COVID-19 Research Project
# Purpose  : Create county level maps 
#############################################

# clear environment 
rm(list = ls())

# load packages 
pkgTest <- function(x) {
  if (!require(x, character.only = TRUE))
  {
    install.packages(x, dep = TRUE)
    if (!require(x, character.only = TRUE))
      stop("Package not found")
  }
}

packages <- c("rgdal" , 
              "ggplot2" ,
              "dplyr" , 
              "foreign" , 
              "maps" , 
              "socviz", 
              "maps", "NLP" ,  
              "usmap", "readxl" , "dplyr" , "viridis")
lapply(packages, pkgTest)

#set working directory 
getwd()
setwd("C:/Users/g1stm01/Dropbox (Research)/stay_at_home_project")

#======================================
# Prepare County Lat/Long Data 
#======================================

# us county map data 
library(usmap)
counties.geo = us_map(regions = "counties")
map_counties = ggplot() + geom_path(data = counties.geo , aes(x = x, y = y , group = group))
map_counties

#======================================
# Prepare State Lat/Long Data 
#======================================

library(usmap)
states.geo = us_map(regions = "states")
map_states = ggplot() + geom_path(data = states.geo , aes(x = x, y = y , group = group))
map_states


#=================================================
# Stay at home data - state level map  
#=================================================

sah_state = read.dta("data/clean/stay_at_home/stay_at_home_stmaps.dta")

remove(states.merged)
states.merged = left_join(states.geo , sah_state, by = "fips" , all.x = TRUE , all.y = TRUE)
sum(is.na(states.merged$st_group))

lbls = vector()
lbls[1] = ""
lbls[1] = (sah_state$st_group_lbl[sah_state$st_group==1])[1]
lbls[2] = (sah_state$st_group_lbl[sah_state$st_group==2])[1]
lbls[3] = (sah_state$st_group_lbl[sah_state$st_group==3])[1]
lbls[4] = (sah_state$st_group_lbl[sah_state$st_group==4])[1]

map_sah_states = ggplot() + 
                 geom_polygon(data = states.merged , aes(x = x, y = y , group = group , fill = as.factor(st_group))) + 
                 geom_path(data = states.geo , aes(x = x, y = y , group = group), color = "black") + 
                 theme_void() + theme(legend.position = "bottom" , legend.text = element_text(size = 14), 
                                      legend.title = element_text(size = 10)) + 
                 guides(fill = guide_legend(nrow = 4)) + 
                 scale_fill_manual(labels = lbls , name = "" , 
                                                    values = c("#08519c" , "#3182bd" , "#9ecae1", "#eff3ff"))

map_sah_states
map_name = paste('maps/state_stay_at_home', '.png' ,  sep = "" )
print(map_sah_states)
dev.copy(png, map_name , width = 1000 , height = 600) 
dev.off()



#=================================================
# Stay at home data - county level map 
#=================================================

# clean stay at home data 
stayathome = read.dta("data/clean/stay_at_home/stay_at_home_maps.dta")

library(dplyr)
stayathome = stayathome %>% rename(fips = county_fips)
colnames(stayathome)
stayathome$fips = as.character(stayathome$fips)

attach(stayathome)
stayathome$len = nchar(fips)
stayathome$zero = ifelse(stayathome$len==4,0,NaN)
stayathome$fips = ifelse(stayathome$len==4,paste(stayathome$zero,fips  , sep = ""), 
                              fips)
detach(stayathome)

# merge with counties data 
remove(counties.geo.merged)
counties.geo.merged = left_join(counties.geo , stayathome, by = "fips" , all.x = TRUE , all.y = TRUE)
sum(is.na(counties.geo.merged$effective_date))


# maap - version1 
map_stayathome_cnty = ggplot() + 
                      geom_polygon(data = counties.geo.merged , aes(x = x, y = y , group = group , fill = as.factor(day_id_groups))) + 
                      geom_path(data = counties.geo.merged, aes(x = x , y = y, group = group) , color = NA ) + 
                      geom_path(data = states.geo , aes(x = x, y = y , group = group), color = "black") + 
                      ggtitle("US counties : Effective Date of Stay At Home Orders") +
                      theme_void() + theme(plot.title = element_text(hjust=0.5)) + 
                      scale_fill_brewer(palette = "OrRd" , direction = -1, name = "Effective Date", 
                                        labels = c("Mar 17 - Mar 19" , "Mar 20 - Mar 22" , "Mar 23 - Mar 25" ,  "Mar 26 - Mar 28", 
                                                   "Mar 29 - Mar 31" , "Apr 1 - Apr 3" , "Apr 4 - Apr 6" , "N/A"))  
map_stayathome_cnty
map_name = paste('maps/county_stay_at_home', '.png' ,  sep = "" )
print(map_stayathome_cnty)
dev.copy(png, map_name , width = 1000 , height = 600) 
dev.off()



# map - version2 
map_stayathome_cnty = ggplot() + 
  geom_polygon(data = counties.geo.merged , aes(x = x, y = y , group = group , fill = day_id), show.legend = TRUE) + 
  geom_path(data = counties.geo.merged, aes(x = x , y = y, group = group) , color = NA ) + 
  geom_path(data = states.geo , aes(x = x, y = y , group = group), color = "black") + 
  theme_void() + theme(plot.title = element_text(hjust=0.5 , size = 18)) + 
  theme(legend.text = element_text(size = 13), 
        legend.title = element_text(size = 15), legend.position =  "bottom" , legend.key.width = unit(2.5,"cm")) + 
  scale_fill_gradientn( colours = c( '#67000d' , '#99000d', '#a50f15', '#de2d26', '#cb181d', 
                                     '#ef3b2c', '#fb6a4a', '#fcbba1', '#fc9272', '#fee0d2', '#fff5f0'),
                        labels = c("Mar 17", "Mar 21" , "Mar 26" , "Mar 31" , "Apr 5" , "Apr 7") ,name = "Effective \n Date" , 
                        limits = c(0, 25) , 
                        na.value = "white")
map_stayathome_cnty
map_name = paste('maps/county_stay_at_home_v2', '.png' ,  sep = "" )
print(map_stayathome_cnty)
dev.copy(png, map_name , width = 1000 , height = 600) 
dev.off()


#=================================================
# Stay at home data - Maps for Gif 
#=================================================

# clean stay at home data 
stayathome = read.dta("data/clean/stay_at_home/stay_at_home_forR.dta")

library(dplyr)
stayathome = stayathome %>% rename(fips = county_fips)
colnames(stayathome)
stayathome$fips = as.character(stayathome$fips)

attach(stayathome)
stayathome$len = nchar(fips)
stayathome$zero = ifelse(stayathome$len==4,0,NaN)
stayathome$fips = ifelse(stayathome$len==4,paste(stayathome$zero,fips  , sep = ""), 
                         fips)
detach(stayathome)

# omit values with NA effective_date 
stayathome = stayathome[!is.na(stayathome$effective_date),]

# get a vector of unique dates in the data 
dates = unique(stayathome$date)
leng_dates = length(dates)

# loop through dates 
k = 1 
while(k <= leng_dates){ 
  
  print(k)
  
  data_small = stayathome[which(stayathome$date == dates[k]),]
  data_small = data_small[which(data_small$stay_at_home == 1),]
  data_small = data_small[,c("fips", 'stay_at_home')]
  
  # merge with counties data 
  merged_fips = left_join(counties.geo , data_small , by = "fips" , all.x = TRUE , all.y = TRUE)

  
 
  map_stayathome_cnty = ggplot() + 
    geom_polygon(data = merged_fips , aes(x = x, y = y , group = group , fill = stay_at_home),  show.legend = FALSE) + 
    geom_path(data = merged_fips , aes(x = x , y = y, group = group) , color = 'NA' ) + 
    geom_path(data = states.geo , aes(x = x, y = y , group = group), color = "black") + 
    ggtitle(dates[k]) + 
    theme_void() + theme(plot.title = element_text(hjust=0.5, size = 25))  + scale_fill_gradient(low = '#a50f15' , high = '#a50f15' , na.value = 'white')
  
 
  map_stayathome_cnty
  map_name = paste('maps/for_gif_date_', k, '.png' ,  sep = "" )
  print(map_stayathome_cnty)
  dev.copy(png, map_name , width = 1000 , height = 600) 
  dev.off()
  
  merged_fips = merged_fips[,-which(names(merged_fips) %in% c("stay_at_home"))]
  k = k + 1
  

} 



#=================================================
# Womply - county level map 
#=================================================

# clean stay at home data 
womply = read.dta("data/clean/womply/womply_formaps.dta")
womply$rev_Other = format(womply$rev_Other , format = "e" , digits = 3)
womply$rev_Restaurant = format(womply$rev_Restaurant , format = "e" , digits = 3)

womply$pc_rev_Other = round(womply$pc_rev_Other, digits = 2 )
womply$pc_rev_Restaurant = round(womply$pc_rev_Restaurant, digits = 2)


# merge with counties data 
rm(merged_womply)

merged_womply = left_join(counties.geo , womply , by = "fips" , all.x = TRUE , all.y = TRUE)

# set functions 
colMax = function(data) sapply(data, max, na.rm = TRUE)
colMin = function(data) sapply(data, min, na.rm = TRUE)

min = min(womply$pc_rev_Other)
max = max(womply$pc_rev_Other)
mid = median(womply$pc_rev_Other)

# create a map - Other 
library(viridis)
map_other = ggplot() + 
            geom_polygon(data = merged_womply , aes(x = x, y = y , group = group , fill = gr_Other),  show.legend = TRUE) + 
            geom_path(data = merged_womply , aes(x = x , y = y, group = group) , color = 'NA' ) + 
            geom_path(data = states.geo , aes(x = x, y = y , group = group), color = "black") + 
            theme_void() + theme(plot.title = element_text(hjust=0.5, size = 25)) + 
          theme(legend.text = element_text(size = 13), 
                  legend.title = element_text(size = 15) , legend.position = "bottom") +
            scale_fill_gradient(low = "white" , high = "chocolate1" , na.value = "white" , breaks = c(0,5,10) , 
                                limits = c(0,10), 
                                labels = c(min , mid , max), 
                                name = "Per-Capita \n Spending ")

map_other
map_name = paste('maps/womply_other', '.png' ,  sep = "" )
print(map_other)
dev.copy(png, map_name , width = 1000 , height = 600) 
dev.off()


# create a map - Restuarant

min = min(womply$pc_rev_Restaurant)
max = max(womply$pc_rev_Restaurant)
mid = median(womply$pc_rev_Restaurant)

library(viridis)
map_restaurant = ggplot() + 
                 geom_polygon(data = merged_womply , aes(x = x, y = y , group = group , fill = gr_Restaurant),  show.legend = TRUE) + 
                 geom_path(data = merged_womply , aes(x = x , y = y, group = group) , color = 'NA' ) + 
                 geom_path(data = states.geo , aes(x = x, y = y , group = group), color = "black") + 
                 theme_void() + theme(plot.title = element_text(hjust=0.5, size = 25)) + 
                 theme(legend.text = element_text(size = 13), 
                       legend.title = element_text(size = 15) , legend.position = "bottom") + 
                  scale_fill_gradient(low = "white" , high = "coral1" , na.value = "white" , breaks = c(0,5,10) , 
                                      limits = c(0,10), 
                                      labels = c(min , mid , max), 
                                      name = "Per-Capita \n Spending") 

map_restaurant
map_name = paste('maps/womply_restaurant', '.png' ,  sep = "" )
print(map_restaurant)
dev.copy(png, map_name , width = 1000 , height = 600) 
dev.off()
